Find spatial configuration of 1ha restoration plots (10,000m2, 100m*100m) at Moore Reef (Cairns). The code takes the Allen Coral atlas geomorphology data (see here for more details), and runs the following process:

  1. extract all Sheltered Reef Slope polygons
set.seed(1)

# data packages
library(tidyverse)
library(units)
library(foreach)
library(future.apply)

# spatial packages
library(sf)

# mapping packages
library(tmap)
library(leaflet)

set_restoration_plot_centres <- function(input = NULL, width = NULL, length = NULL) {
  # Calculate the coordinates of the rectangular polygon
  x <- sf::st_coordinates(input)[1, 1]
  y <- sf::st_coordinates(input)[1, 2]

  # set parameters
  x_min <- x - (width / 2)
  x_max <- x + (width / 2)
  y_min <- y - (length / 2)
  y_max <- y + (length / 2)

  polygon <- sf::st_polygon(list(rbind(c(x_min, y_min), c(x_min, y_max), c(x_max, y_max), c(x_max, y_min), c(x_min, y_min)))) |>
    sf::st_sfc(crs = 20353)

  return(polygon)
}

habitats <- c("Plateau", "Back Reef Slope", "Reef Slope", "Sheltered Reef Slope", "Inner Reef Flat", "Outer Reef Flat", "Reef Crest")
raw_geomorphic <- st_read("/Users/rof011/GBR-coral-restoration/data/Moore-Reef-20230906030257/Geomorphic-Map/geomorphic.geojson", quiet=TRUE) %>%
  sf::st_transform(crs = sf::st_crs("EPSG:20353")) %>%
  dplyr::group_by(class) %>%
  dplyr::mutate(habitat_id = paste0(
    gsub(" ", "_", class), "_",
    sprintf(paste0("%0", ceiling(log10(max(1:length(class)))), "d"), 1:length(class))
  )) %>%
  sf::st_make_valid() %>%
  dplyr::group_by(habitat_id, class) %>%
  summarise(geometry = st_union(geometry)) %>%
  mutate(area = round(st_area(geometry))) %>%
  filter(class %in% habitats) %>%
  drop_na(class)

moore_geomorphic <- st_read("/Users/rof011/GBR-coral-restoration/data/Moore-Reef-20230906030257/Geomorphic-Map/geomorphic.geojson", quiet=TRUE) %>%
  sf::st_transform(crs = sf::st_crs("EPSG:20353")) %>%
  dplyr::group_by(class) %>%
  dplyr::mutate(habitat_id = paste0(
    gsub(" ", "_", class), "_",
    sprintf(paste0("%0", ceiling(log10(max(1:length(class)))), "d"), 1:length(class))
  )) %>%
  sf::st_make_valid() %>%
  #  mutate(habitat_id=as.factor(habitat_id)) %>%
  dplyr::group_by(habitat_id, class) %>%
  summarise(geometry = st_union(geometry)) %>%
  mutate(area = round(st_area(geometry))) %>%
  filter(class %in% habitats) %>%
  drop_na(class)
  1. extract Outer Reef Flat polygons bordering Sheltered Reef Slope polygons (removes exposed Outer Reef Flats and Inner Reef Flats that are windward of Sheltered Reef Slopes) using sf::st_touches() and purrrr::map2().
# find restorable habitats
moore_sheltered_sites <- filter(moore_geomorphic, class %in% c("Sheltered Reef Slope", "Outer Reef Flat")) 
# identify neighbouring pixels
neighbors <- st_touches(moore_geomorphic)

# match pairs
border_rows <- purrr::map2(neighbors, moore_geomorphic$class, ~ {
  any(moore_geomorphic$class[.x] == "Sheltered Reef Slope") & (.y == "Outer Reef Flat")
})

# Filter rows
indices <- which(unlist(border_rows))
result <- moore_geomorphic[indices, ]

# combine outer and sheltered
moore_restorable <- rbind(
  filter(moore_geomorphic, class %in% c("Sheltered Reef Slope")),
  result
)

# make a single polygon set for seeding with points. 
# merge classes, union, filter by >15,000m2
moore_restorable_pooled <- moore_restorable %>% 
  select(-class, -habitat_id) %>%
  summarise(geometry = st_union(geometry)) %>%
  st_cast("POLYGON") %>%
  mutate(area = round(st_area(geometry))) %>%
  filter(area > set_units(15000, m2))
  1. add a negative buffer (-50m) to polygons using sf::st_buffer, drop null geometries (by assuming that centroids of 1ha plots are not <50m from reef boundaries the buffer speeds up placement of random 1ha polygons)
moore_buffered <- moore_restorable_pooled %>% 
  st_buffer(set_units(-50, m)) %>%
  filter(!st_is_empty(.))
  1. seed remaining polygons with 100k random points, apply using function and multicore to speed up the process:
  1. expand a 100 x 100m polygon around each point
  2. retain if polygon overlaps 100% with underlying reef habitat using sf::st_within(), remove if overlap
  3. iterate for n points using future.lapply (~600 seconds for 100k polygons)
set_plot <- function(input = NULL, width = NULL, length = NULL, center = NULL) {
# Calculate the coordinates of the rectangular polygon
x <- sf::st_coordinates(input)[1, 1]
y <- sf::st_coordinates(input)[1, 2]
# set parameters
x_min <- x - (width / 2)
x_max <- x + (width / 2)
y_min <- y - (length / 2)
y_max <- y + (length / 2)
polygon <- sf::st_polygon(list(rbind(c(x_min, y_min), c(x_min, y_max), c(x_max, y_max), c(x_max, y_min), c(x_min, y_min)))) |>
sf::st_sfc(crs = 20353)
return(polygon)
}

moore_buffered_large <- filter(moore_buffered, area>units::set_units(15000, m2))

generate_random_polygon <- function(input) {
  random_point <- st_sample(moore_buffered_large, 1)
  random_polygon <- set_plot(random_point, 100, 100)
  check_overlap <- st_within(random_polygon, moore_restorable_pooled)
  if (isTRUE(check_overlap[1] >= 1)) {
    return(random_polygon)
  }
}

num_iterations <- 100000 # set iterations per polygon
valid_polygons <- list() # create list

# Use future_lapply to parallelize the iterations
plan(multisession, workers = 11)
valid_polygons_moore_buffered <- future_lapply(1:num_iterations, future.seed = 1, function(i) {
  generate_random_polygon() # filter lower boundary
})

# Unlist the results to get a list of valid polygons
valid_polygons_moore_buffered_final <- unlist(valid_polygons_moore_buffered, recursive = FALSE) |> st_sfc(crs=20353)

valid_polygons_moore_buffered_final_stc <- valid_polygons_moore_buffered_final %>% 
  st_as_sf(crs=20353) %>% 
  mutate(polygon.ID = seq(1:length (valid_polygons_moore_buffered_final)))

Map restoration plots

The white plot outlines are a small subset of potential polygon outlines (200 from 100k), pink plots represent manually selected non-overlapping plots.

Select layers using the layercontrol on the left, use [ ] for full-screen viewing.

library(tmap)


# for now manually subset polygons
subsetpolygons <- c(

c(35, 59, 188, 186),
c(142, 13, 190),
c(144,174),
c(54, 169, 66, 140, 60, 9)
)

valid_polygons_moore_buffered_final_nonoverlap <- valid_polygons_moore_buffered_final_stc |> filter(polygon.ID %in% subsetpolygons)

# visualise with thematic maps
map <- tm_view() +
#tm_tiles("Esri.WorldImagery", group = "[Satellite map]", alpha = 0.5) +

#---------- all habitats--------@
tm_shape(raw_geomorphic |> mutate(class=as.factor(class)),  name = "Reef habitats") +
  tm_borders(col = "black", lwd = 0.2) +
  tm_fill("class", name = "Benthic habitats", Title = "Benthic habitats", id="area", palette = c("Plateau" = "lightskyblue1", "Back Reef Slope" = "darkcyan",
                                                                    "Reef Slope" = "darkseagreen4", "Sheltered Reef Slope" = "darkslategrey",
                                                                    "Inner Reef Flat" = "darkgoldenrod4", "Outer Reef Flat" = "darkgoldenrod2",
                                                                    "Reef Crest" = "coral3"), alpha = 0.6) +
#---------- restorable habitat--------@
tm_shape(moore_restorable |> mutate(class=as.factor(class)), id="area", name = "Restorable habitats") +
  tm_borders(col = "white", lwd = 0) +
  tm_fill("class", title="Restorable habitats", name = "Restorable habitats", id="area", palette = c("Sheltered Reef Slope" = "darkslategrey",  "Outer Reef Flat" = "darkgoldenrod3"), alpha = 0.6) +

#---------- buffered habitat--------@

tm_shape(moore_buffered, name = "Buffered habitats") +
  tm_borders(col="darkred", lwd=0.2) +

#---------- restoration plots--------@

 tm_shape(valid_polygons_moore_buffered_final[1:200], name = "Restorable plots") +
   tm_borders(col = "white", lwd = 0.5, alpha=0.2,) +
   tm_fill(col = "white", alpha = 0, name = "Polygon ID", id="polygon.ID") +
  
#---------- restoration plots--------@
   tm_shape(valid_polygons_moore_buffered_final_nonoverlap, name="Non-overlapping plots", is.master=TRUE) +
   tm_fill(col = "deeppink3", alpha = 0.1, name = "subset.polygons", id="subset.polygons") +
   tm_borders(col="deeppink3", lwd=1) +
 
  
#---------- scale--------@
tm_scale_bar(width=200)

  map |> tmap::tmap_leaflet() |>
    leaflet::addProviderTiles('Esri.WorldImagery', group = "Satellite map", options=leaflet::providerTileOptions(opacity=0.6, maxNativeZoom=18,maxZoom=100)) |>
#    leaflet::addProviderTiles('Esri.WorldTopoMap',  group = "<b> [Restoration]</b> habitats", options=leaflet::providerTileOptions(maxNativeZoom=19,maxZoom=10)) |>
    leaflet::addLayersControl(position="topleft", overlayGroups=c(
      "Satellite map", 
      "Reef habitats",
      "Restorable habitats",
      "Non-overlapping plots",
      "Buffered habitats",
      "Restoration plots (1ha)"
      ),
                              options=leaflet::layersControlOptions(collapsed = FALSE)) |>
     leaflet::hideGroup(c("Reef habitats")) |> 
    leaflet.extras::addFullscreenControl(position = "topleft", pseudoFullscreen = FALSE)

Summary statistics

Habitat Area TwentyPlots
Total coral area 9637437 2.1%
Total restorable area 2686748 7.4%